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Chest physiotherapy is an empirical technique used to help secretions to get out of the lung 
whenever stagnation occurs. Although commonly used, little is known about the inner mechanisms 
of chest physiotherapy and controversies about its use are coming out regularly. Thus, a scientihc 
validation of chest physiotherapy is needed to evaluate its effects on secretions. 

We setup a quasi-static numerical model of chest physiotherapy based on thorax and lung physiol¬ 
ogy and on their respective biophysics. We modeled the lung with an idealized deformable symmetric 
bifurcating tree. Bronchi and their inner fluids mechanics are assumed axisymmetric. Static data 
from the literature is used to build a model for the lung’s mechanics. Secretions motion is the conse¬ 
quence of the shear constraints apply by the air flow. The input of the model is the pressure on the 
chest wall at each time, and the output is the bronchi geometry and air and secretions properties. 

In the limit of our model, we mimicked manual and mechanical chest physiotherapy techniques. 

We show that for secretions to move, air flow has to be high enough to overcome secretion resistance 
to motion. Moreover, the higher the pressure or the quicker it is applied, the higher is the air flow 
and thus the mobilization of secretions. However, pressures too high are efficient up to a point 
where airways compressions prevents air flow to increases any further. Generally, the hrst effects of 
manipulations is a decrease of the airway tree hydrodynamic resistance, thus improving ventilation 
even if secretions do not get out of the lungs. Also, some secretions might be pushed deeper into 
the lungs; this effect is stronger for high pressures and for mechanical chest physiotherapy. Finally, 
we propose and tested two adimensional numbers that depend on lung properties and that allow to 
measure the efficiency and comfort of a manipulation. 
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I. INTRODUCTION 


Lung secretions are used to protect the lungs against pollutants or infections, by capturing airborne particles or cells 
that may enter the lungs. Secretions form an heterogenous gel-like fluid consisting mostly in water and biopolymers 
M- Secretions are motioned upwards thanks to two natural mechanisms. First, the cilia located on the bronchi 
walls, whose motion pushes the secretions upward the bronchial tree, this phenomenon is called mucociliary clearance 
[iniiisiEQ]). The second mechanism is cough, it involves high air flows into the bronchi [3l [311 El] ■ These mechanisms 
can however fail. For example, cough or cilia motion are not mature in young children, they can also been altered 
by pathologies or age. As a consequence, when secretions are too thick or overproduced, they can become stalled 
and bronchial infections may occur. Chest physiotherapy is used to treat patients whose secretions natural clearing 
mechanisms are altered. The techniques depend on the pathology and on the patient characteristics, typically her/his 
age. Bronchiolitis, asthma, chronic obstructive pulmonary disease or cystic fibrosis are typical diseases treated with 
dedicated chest physiotherapy techniques. 

Chest physiotherapy is based on external maneuvers on the chest. It relies on a basic concept: the pressures 
applied by the practitioner on the chest transmit into the lungs and deforms the bronchi. Bronchi’s volume changes 
create an air flow inside the lung and the mechanical constraints applied on the secretions by the air flow affect the 
secretions and have the potential to induce motion. Secretions are gel-like materials which are motioned only if they 
are submitted to a stress sufficiently strong. In manual chest physiotherapy, flow and lung volume are modulated by 
the practitioners during the manipulations in response to the feedback they feel from the patient. Mechanical devices 
based on this concept have also been developed [i[iaii7|. This concept represents however a very simplified view of 
the lung and thorax biomechanics and chest physiotherapy remains as of today a discipline which is mostly empirical. 
Maneuvers are efficient but the inner mechanisms involved in their efficiency remain either unknown or non proved 
in a scientific way. Many chest physiotherapy techniques exist, and their definition and whether they should be used 
or not are decided during country wide consensus meetings. Consequently, the available pool of maneuvers can differ 
from one country to the other, mostly because their efficiency are difficult to compare [I^ 128) . 

Historically, most of the maneuvers were based on applying relatively strong pressures on the chest, bringing 
discomfort and, if ill performed, to risks to the patient, most particularly to baby and children. Thus, the use 
of manual chest physiotherapy for baby suffering bronchiolitis - a well known and wide spread disease - has been 
questioned because of the possible risks induced by the manipulations HHTj. Recently, new techniques focusing on 
patient comfort and using low pressures on the chest, low air flows in the lungs and induced cough have emerged and 
have spread quickly in physiotherapists’ community. The efficiency of these new low flow techniques is however as of 
today still controversial, because chest physiotherapy is based on the hypothesis that secretions can be mobilized (i.e. 
put to motion) only with high air flows. 

Chest physiotherapy has reached a state where it needs to be supported by a strong scientific background. The 
different questions raised by physiotherapists would benefit from a better understanding of the coupled biomechanics 
of chest and lungs. An integrated model of chest physiotherapy that is able to uncover the effects of hand pressures 
on secretions mobilization would help the practitioners to chose and define knowledgeably the techniques they are 
using. Most particularly, the role of the bronchial tree geometry play an important role on the air fluid dynamics 
[smillllllig and may have an important role on secretions mobilization. Hence, in a previous paper [20) . we 
studied the role of the interaction between air flow, secretion thickness and bronchial tree geometry assuming it was 
rigid. This first theoretical frame contained some of the core elements involved in chest physiotherapy and we showed 
that air flow needed to be high enough for secretions to be motioned and that applying the same given air flow in the 
lung leads to a stationary distribution of mucus, i.e. a constant flow manipulation might see its efficiency decrease 
over time. These first results indicate that air flow needs to be modulated towards higher air flows in order to be able 
to mobilize mucus all along the manipulation. Another method for increasing the air flow induced stress might be to 
decrease the volume of the bronchi while keeping the air flow rate constant. This is the idea used, for example, in low 
volume manipulations [28j or in the mechanical chest physiotherapy technique of chest compression IldTj. In order 
for our model to be able to capture such behaviors and to give some qualitative and comparative predictions, it was 
necessary to add to our previous model how lung, and its subunits bronchi and acini, deform during the manipulations. 
Moreover, integrating lung deformations would improve globally the quality of the predictions of our model. 

Thus, in this work, we propose a new model that integrates lung deformations. The model input is an homogeneous 
pressure applied on the chest, and the model output is the air flow, the secretions motion and distribution and the 
lung configuration (airways sizes, tissue pressure, etc.). In section |Hj we describe the hypotheses of the model and in 
section m we present and discuss the predictions of the model in the case of idealized manipulations of manual chest 
physiotherapy and of different devices of mechanical chest physiotherapy. 
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II. METHODS 

This study is based on three separate sub-models: one for the lung, one for the thorax and one to mimic chest 
physiotherapy manipulations. We build a minimal models which includes only the core properties we identified as 
playing an important role in chest physiotherapy. This choice was mainly made in order to be able to interpret correctly 
the interactions between these core phenomena, it was also a way to limit the computation times of the numerical 
studies. As a consequence, our model can perform qualitative and comparative predictions, but not quantitative 
predictions. 

This work assumes that lung tissue mechanics and lung generations are homogeneous. These choices were made 
to keep the model tractable in term of physical behavior. Thus we chose to neglect phenomena such as hydrostatic 
pressures, inhomogeneous pressure distribution in the lung tissue or differences between airways or alveoli belonging 
to the same generations. The details of the model are given in the following subsections. 


A. Terminology. 

Below is a list of definitions for some terms or variables used in this work. 

FRC: functional residual volume, the volume of the lung at the end of a normal expiration. 

Generation: the generation of an airway is an index accounting for the number of bifurcations in the path from the 

tree root (trachea) to that airway (tree root is zero, deepest alveoli are twenty-three). We denote z{i) the 

generation index of airway i. 

Lung volume: lung volume is denoted Vl it is the sum of the tracheobronchial tree volume Vtbt and of the acini volume 

Kc- 

Tree structure: for an airway indexed i, its parent branch is indexed p, and its daughter branches are indexed di and 
d 2 - Although p, di and ^2 are function of z, we do not make i appear for the sake of notations simplification. 

1 cmH20: 1 cmH20 is equivalent to 98.07 Pa. 

iz) 

•. air flow in an airway in generation 2 . 

(z) 

/• aiucus flow in an airway in generation z. 

(z) 

alveolus connected to an airway in generation z. 

G: length of the bronchi from generation index z. 

max(a, b): maximum function, max(a, b) = a is a > b and max(a, b) = b otherwise. 

i Z^ 

number of alveoli connected to an airway in generation z (0 for airways which are not alveolar ducts, 58 for 
alveolar ducts, i.e. whose generation is deeper than 17). 

Pext- external chest wall pressure. It is assumed constant everywhere on the chest. 

(2) 

P^i^'. air pressure at mid-airway length for airway i. 

PtissueiVL)- mean lung tissue pressure when the volume of the lung is Vl- This data comes from static lung volume- 
pressure relationships measured by Agostoni et al. [T]. 

r^a'^: air lumen area radius. 

airway lumen area. 

i Z^ 

Sa ■ air lumen area of bronchus i. Air lumen area is assumed disk shaped. 

(z) 

: airway i lumen area. Airway lumen area is assumed disk shaped. 

VL.rsi^P)- Lung volume when it is submitted to a pressure difference AP between the thorax and the bronchi 
(same pressure into the whole lung, no air movement). This quantity is built from static lung volume-pressure 
relationships measured by Agostoni et al. [T]. 
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z: generation index in the tree. 

Note that the subscript {z) corresponds to the airway generation associated to the data. In the case of properties 
that are applicable to any airway, we may drop this subscript for the sake of formulas simplification. 


B. Lung’s geometry 


The geometrical model of the lung used in this work is scaled to have a total lung capacity (TLC) of 6.5 L, a vital 
capacity (VC) of 5 L, a residual volume (RV) of 1.5 L and a tidal volume at rest of 0.5 L [3132]. According to [3, 
functional residual volume (FRC) is reached at 35 % of VC, thus FRC = RV + 0.35 x VC is 3.25 L. 


1. Bronchial tree model. 


The bronchial tree is modeled as an idealized airway tree with twenty-three generations. The generation index starts 
from generation zero (trachea) to generation twenty-two (deepest alveolar duct) and is incremented by one at each 
bifurcation. Although bronchial bifurcations are slightly asymmetrical 13111 muni, we approximate the bronchial tree 
with a symmetric airway tree for our pioneering model to remain tractable. We assume that each airway bifurcates 
into two smaller identical airways, implying that all the airways belonging to a same generation are identical. 

The sizes and lengths of the conducting airways (seventeen first generations) are chosen from Lambert’s data [16) . 
the different parameters are recalled in table |l| The airways are modeled as cylindrical tubes. We assume that mucus 
distribution in an airway is a constant thickness layer along the cylinder wall and that mucus cannot reach the alveoli. 


z 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

1 (cm) 

12.00 

4.76 

1.90 

0.76 

1.27 

1.07 

0.90 

0.76 

0.64 

0.54 

0.47 

0.39 

0.33 

0.27 

0.23 

0.20 

0.17 

d (cm) 

1.671 

1.1815 

0.8735 

0.670 

0.525 

0.399 

0.3095 

0.241 

0.192 

0.151 

0.119 

0.096 

0.080 

0.070 

0.0615 

0.055 

0.0495 


Table I. Lung’s data for the geometry of the conductive airways at FRC = 1.75 L (end expiration volume), from Lambert et 
al [16]: z is the generation index; I is the airway length; d is the airway diameter. 


2. Alveoli and acini geometrical models. 

Alveoli are assumed to be spherical and to stem from airways walls from the seventeenth generation. The total 
number of alveoli Na is documented in the literature [33] and we used Na = 480 10®. FRC volume can be decomposed 
as the volume of the bronchial tree plus the volume of the alveoli. Our model computations give a FRC volume for 
the bronchial tree of about 0.67 L and a FRC total alveoli volume of about 2.58 L. Consequently, the volume of an 
alveolus at FRC used in our model is Va = 5.37 10® nm? (diameter of 217 ixm). This value is slightly higher than the 
value measured in |26j (4.2 10® /rm®, i.e. a diameter of 200 /rm), the discrepancy may arise from the measures being 
done in a lung’s volume that is not FRC. The exchange surface in the lungs consists in the total alveoli surface and 
can be estimated in our model for FRC to ^ 70 m^. 

Alveoli are stemming from bronchi walls in the acini units, that corresponds to the last sixth generations of the 
lungs. In our model, each acinus has 2® — 1 = 63 alveolar ducts. On the contrary of bronchi in the seventeen first 
generations whose diameters and lengths decrease when generation index increases, the bronchi in the acini, called 
alveolar ducts, keep roughly similar size whatever their generation [32) . Thus, we assume that the number of alveoli 
is homogeneously distributed on the alveolar ducts walls, and to reach a total of 480 10® alveoli, the amount of alveoli 
per alveolar ducts in the acini is 58, a bit higher than the range from 20 to 40 reported in [33]. It remains however 
fully in the range of physiological patterns, since there is quite large variation between individuals, as reported for 
example in [26) . 


C. Lung’s mechanics 


Our goal is to be able to compute the volume of the airways and alveoli as a function of two input pressures: the 
pressure applied on the thorax Pext and the fluid pressures inside the airways Pk,: {z is the generation index). In 
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our model, the pressure applied on the thorax may either come from muscles action, like in ventilation, or from chest 
physiotherapy manipulation; the pressure inside the airways is a consequence of the fluid mechanics in the tree. For 
the sake of simplification, we assume that the pressure in the lung tissue Pussue to be homogeneous throughout the 
lung. We assume that lung tissue pressure Pussue depends on the lung volume only and is computed from the static 
volume-pressure relationship from Agostoni et al work mm , see hgure (dashed curve). Similarly, we assume that 
bronchi and alveoli mechanics do not depend on their location in the lung. 

To model the lung mechanics, we divided the tree in two parts: the conductive airways and the acini. Mechanics of 
airways is based on the transmural pressure-airways section relationships proposed by Lambert et al |16j . Mechanics 
of the acini is based on the total respiratory system lung volume-pressure relationship from Agostoni et al mm , see 
figure (continuous curve). 



Figure 1. Data (static) from Agostoni et al. [I]. Continuous curve: total respiratory lung volume-pressure relationship, used 
in our work to model acini mechanics. Dashed curve: Lung tissue pressure, used in our work to compute bronchi transmural 
pressure. 


1. Conductive airways mechanics (up to the 17th generation). 


To compute airways diameters, we used data from Lambert et al. m- Lambert et al built data based static 
relationships of bronchi’s airways lumen area versus their transmural pressure for the seventeenth first generations. 
Their estimations are given per generation. For an airway belonging to generation z with a transmural pressure AP, 
they proposed that the airway lumen area ^'^(Ap) is a sigmoid function of transmural pressure: 


5.(AP) = 


( \ (^) 

^miz) whenAP<0 

(^1 - (1 - ao(-2))(l - A^iz) when AP > 0 


( 1 ) 


with Pi(z) = and P 2 {z) = —n 2 {z)^-^^j^■ Am represents the maximal possible lumen area. The quantity 

a^Am is the airway area when transmural pressure is 0. At functional residual regime (FRC, end of normal expiration), 
the transmural pressure is about 500 Pa. The values used for the different parameters in the previous formulas are 
detailed in table |llj the dependence of the bronchi surface area with transmural pressure are plotted on figure 


z 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

Oio 

0.882 

0.882 

0.686 

0.546 

0.450 

0.370 

0.310 

0.255 

0.213 

0.184 

0.153 

0.125 

0.100 

0.075 

0.057 

0.045 

0.039 

Oq (xlO”'” Pa~^) 

1.1 

1.1 

5.1 

8.0 

10.0 

12.5 

14.2 

15.9 

17.4 

18.4 

19.4 

20.6 

21.8 

22.6 

23.3 

23.9 

24.3 

ni 

0.5 

0.5 

0.6 

0.6 

0.7 

0.8 

0.9 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

n2 

10.00 

10.00 

10.00 

10.00 

10.00 

10.00 

10.00 

10.00 

10.00 

10.00 

10.00 

9.00 

8.00 

8.00 

8.00 

7.00 

7.00 

Am (cm’^) 

2.37 

2.37 

2.80 

3.50 

4.50 

5.30 

6.50 

8.00 

10.20 

12.70 

15.94 

20.70 

28.80 

44.50 

69.40 

113.00 

180.00 


Table II. Lung’s data for the static mechanical behavior of the conductive airways, from Lambert et al [16]. z is the generation 
index, see equation Q for the definitions of the other numbers oq, Oq, ni, n 2 and Am- 
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Figure 2. Dependence of bronchi surface area (normalized) with bronchi transmural pressure. The numbers correspond to the 
generation index and the red dashed line represents bronchi state at FRC. 


Airways transmural pressure is the difference between the pressure in the lung’s tissues around the bronchi and the 
bronchi inner pressure due to fluid mechanics in the lung. The fluid mechanics in the airway tree is detailed below. 
We assume that tissue pressure depends on the volume of the lungs only; it is computed from the volume - lung 
pressure static relationship measured by Agostoni et al [T]. Using our terminology, the lumen area of a bronchi 
in generation z {z < 16) is 


= S, (PussueiVL) - ( 2 ) 

We recall that PussueiVL,) is given by the volume-pressure relationship from Agostoni et al. [T], see figure(dashed 
curve). To simplify, we will assume that Pj^^l is evaluated in the middle of the bronchus z. Moreover, the pressure 
- lumen area relationships assumes that airway walls are always at equilibrium with their corresponding transmural 
pressure, thus bronchi mechanics in our model is quasi-static. From these laws, we can derive the volume of the 
tracheobronchial tree, 


Vm{VL,Pa^r) = £ 2 ^ 5 , (PussueiVL) - P^jl) 
2=0 


(3) 


2. Alveolar ducts and alveoli mechanics (last six generations). 

The last sixth generations of the lung correspond to the acini. Each acinus in our model consists in 2® — 1 = 63 
alveolar ducts and each alveolar duct is connected to 58 alveoli (see model geometry section). An alveolar duct with 
its corresponding 58 alveoli will be referred to as an alveolar duct unit. There is 2^^ x 63 alveolar duct units in our 
model. Lambert et al data does not cover those generations, thus we used another approach to approximate the 
mechanical behavior of an alveolar duct unit. 

The pressure difference between the inside of an alveolar duct unit in generation z and the body surface is P))i(.—Pext- 
Consequently, the volume the lung would have under static condition with the same pressure difference in the whole 
tree is given by the respiratory system pressure - volume relationship from [T], see continuous curve on figure ^ 
However, the lung volume measured by Agostoni et al includes the tracheobronchial tree volume, thus to reach the 
volume of acini only, it is necessary to subtract from the lung volume the volume of the tracheobronchial tree as it 
is in the configuration of Agostoni et al measurements: static and the glottis open. As a consequence, the reference 
pressure (= atmospheric pressure) spans the whole tree. Hence, the volume of an alveolar duct unit in generation 
index z in the tree is: 


'^adu 


( pC^) p A 

air ’ ^ e.xt J 


1 


217 X 63 


(^L,rs{Pail Pext) Vtbt (Vl ,r s {Pail Pext),0^'^ 


(4) 


Now, to compute alveolar ducts and alveoli volume, we assume their respective volume ratio in an alveolar duct 
unit remains unchanged whatever the lung volume. We assume the length of alveolar ducts to be Lad = 0.7 mm El- 







7 


The quantity a represents the volume ratio of an alveolar duct into an alveolar duct unit. Hence, for airways in a 
generation z in the acini (z > 16): 

= aVadu /Lad (alveolar duct lumen area) 

’vify = a)vadu(^Pa/ryPext^ (alveolus volume) ^ 

We compute a from FRC state and assume it is a constant: a = 0.17. Finally, we can compute the total volume of 
the acini in the tree, 

22 

Vac{Pa^r,Pext) = {Patl Pext) ( 6 ) 

z=17 


3. Computation of the lung volume Vl ■ 

To compute the conductive airways lumen areas with equation we need to compute the volume of the lung Vl ■ 
Acini volume is known through equation Q and lung volume is the sum of acini volume and of conductive airways 
volume. Consequently, lung volume Vl depends on Pair and Pext and is the solution of the equation: 


Vl = Vac{Pair : Pext ) + VtM{VL ; Pair ) (7) 

Once computed, lung volume Vl is used to compute the tracheobronchial tree lumen areas using equation (§. 
Consequently, the airway lumen area of generations z < 16 expressed in equation (§ can be reformulated as a 

function of Pext and Pair- (^Pussue (VL(Pext, Pa^r)) “ Pj/r) ■ Merging with the equation on the first line of 

equations ([^, one can summarize the airway mechanics into a vectorial form: 

S = H{Pext, Pair) (8) 


with 


H{PexuPa^r)^^'^ = 


{^tissue {VL{PexuPa^r)) - for 0 < Z < 16 

OiVadu (^Pa/l,Pext^ /Lad for 17 < Z < 22 


D. Air and mucus mechanics. 

In the previous section, we defined a model that mimics how airways and alveoli are behaving when they are 
submitted to an external pressure Pext and an internal fluid pressure Pair- Whereas the external pressure Pext is 
an input of the model, fluid pressures in the airways Pair are outputs. Fluid pressures and flows come from the 
interactions between airways deformations and fluid mechanics of air and mucus. 

We will distinguish two lumen areas: the bronchi lumen area, whose surface is measured by the quantity S/ ^ in 

(z) 

generation z, and the air lumen area, whose surface is measured by the quantity Sa in generation z. The difference 

(z) (z) 

S/ — Sa represents the surface of the airway section filled with mucus. 


1. Air mechanics. 


(z) 

In this section, we will derive the equations on the air lumen area surface Sd - Air is assumed incompressible and 
Newtonian. Its density is pa = 1 kg-m~^ and its viscosity is Ha = 1-8 10“® Pas. We assume that air motion in a 
bronchus is at low regime with fully developed profiles as in |20| . The rate of air volume change in an airway is equal 
to the air flow getting in the airway - or getting out the airway, in which case we consider there is a negative airflow 
getting in the airway. Air flow getting in an airway in generation z comes from its two daughter branches standing in 
generation z + 1. The conservation law for air in generation z is: 


H - 



— n 


V) 

alv 


$ 


V) 

alv 


dt 


(9) 



The minus sign means that air flow coming in the tree is modeled as negative. Note that we could also have computed 
the balance with the air flow coming in from the parent branch instead of the daughter branch, however this would 
lead to the similar equations once rewritten at the tree level. The number corresponds to the number of alveoli 
connected to an airway in generation z, it is equal to 0 if the airway belongs to generations zero to sixteen and is equal 
to 58 if it belongs to a deeper generation, see geometry section upwards. corresponds to the part of the airway 

(z) (z) 

air flow that is induced by one alveolus belonging to generation z. is equal to the rate at which the volume 
changes in time. It is computed from equation 


dv 


(^) 

alv 


dt 


= 


(^) 

alv 


( 10 ) 


Because of our hypotheses, the pressure drop per unit length in an airway, called C = ^ depends on the generation 
index only. in an airway of generation z can be computed from the air flow in the bronchus and from air 
and mucus lumen areas, respectively denoted Sa^'^ and 

( 11 ) 

F computation is based on low regime and fully developed fluid dynamics hypotheses, as in our previous work [20) . 
Details about how F is computed are recalled in Appendix]^ Since pressure reference is that of atmospheric pressure, 
the pressure at the root of the tree (trachea) is zero. The pressure drop in an airway of generation z is and 

the pressure in generation z is the sum of the pressure drops of the airways connecting the tree root (trachea) 
to a generation z airway: 


p(^) _ (^(2) 

^air ^ 



g=o 


( 12 ) 


(z) 

The 1/2 in the first term means that 


is air pressure at mid-bronchus length. 


2. Mucus mechanics. 

As in [20], we assume that secretions are a Bingham fluid: it is solid if its inner shear stress a is lower than a 
stress threshold (Tq and liquid with viscosity otherwise. In this work, we assume uq = 0.1 Pa and pm = 0.1 Pa.s, 
see [20] • Secretions density is that of water, i.e. 1000 kg/m^. Mucus is motioned by air flow when the shear forces 
induced by air flow on mucus exceed the shear threshold cto- Mucus can have three states: either fully solid, fully 
liquid or both solid/liquid. The state of mucus depends on a radius ro = |2cro/C'| where C is the pressure drop per 
unit length in the branch [20]. Mucus is solid if > rj,, liquid if tq < ra- If < tq < r^, then mucus is solid between 
radius and Tq and liquid between radius tq and Vf,. 

Mucus volume change in an airway in generation z is the balance between mucus input / and output ‘h/j p. 
Mucus volume in a generation z airway is the difference between airway total volume and the volume occupied by 
air: L/ — Sa Tj, • Since the airway length Lb is assumed constant, the time derivative of this equation leads to 

the conservation law for mucus: 


dS, 




dt 





- 1 $ 


m,0 I 


(13) 


(z) 

The flow of mucus q that gets out of an airway through its interaction with air flow can be computed from air 

(z) 

properties in the airway. Mucus flow q in an airway in generation z is a function of the pressure drop per unit 
length in the airway and of the bronchi and air lumen areas of the airway: 

= (14) 

Details on how to compute G are given in Appendix the method is the same as in [20] ■ Mucus flow in airways 
occurs only if the local shear stress induced by air flow on mucus overcomes the mucus yield stress (Tq. The flow of 
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mucus that gets in a generation z airway comes either from the two daughter airways in generation z + 1 or from the 
parent airway in generation z — 1, depending on mucus flow direction. This is reflected by the sign of mucus flow: 

^ max (o, - 2min (o, (15) 

If mucus flows in the two daughter airways (generation z +1) are negative, then mucus is going up the tree and 

thus accounts twice for mucus inflow in an airway of generation z. If mucus in the parent airway (generation z — 1) 
is positive, then mucus in that airway is going down and accounts half for mucus inflow in an airway of generation 
z. If the flow in one of these neighboring airways are of the opposite sign as stated, then it does not account for the 
mucus inflow in an airway of generation z. 


E. Solving equations. 


The model consists in the set of equations (|8]), ([^, ([^, ([^, ([^, ([^, ([^, ([^. The system of equations is 
implemented numerically in C++ using a vectorial formulation based on Eigen 3.2.1 (http://eigen.tuxfamily.org/) 
and OpenMP (http://openmp.org/) multiprocessing. It is solved thanks to a full implicit first order method for 
differential equations ([^, ([l0| and (13). A Newton method is used to solve equation ^ in order to compute the 
function H, see equation (|^. Details of the mathematical and numerical implementations are given in Appendix 
In order to be able to capture the effects of a 20 i/z oscillating external pressure, the typical time step in our 
simulations is 5.10“^ s. The time needed to compute one time step on 6 cores Xeon 5645 ranges from 0.1 to 0.2 s. 
In this study, the typical manipulation mimicked by our model covers about 200 s, with a total computation time of 
4000 to 8000 seconds. 


III. RESULTS/DISCUSSION 


The model we describe in the previous sections allows to test how the changes in airway diameters affect air and 
mucus distributions in the tree. The model contains the minimal physics that we hypothesize to be main drivers 
of chest physiotherapy: the tree structure, linear air fluid mechanics, Bingham fluid mechanics to model secretions 
behavior and quasi-static deformable airways. In this section, we present the model predictions with pressure inputs 
mimicking chest physiotherapy. In section [ill A[ we define the initial state which is used in all the next simulations. 
This initial state is not representative of any specific pathology, it only reflects a lung with an excess of secretions near 
the eighth generation. In section IIIB[ we present how our model is responding to inputs that mimic manual chest 
physiotherapy. Next, in section C| we present how our model is responding to inputs that mimic high frequency 
chest wall oscillations. We study most particularly the role of an added static pressure on the chest. In order to 
be able to compare both the efficiency and comfort of these manipulations we define in sections |IIID| and |III E| two 
numbers that quantify manipulations efficiency and comfort. 


A. Mimicking chest physiotherapy manipulations 


We will test different input pressures mimicking manual or mechanical chest physiotherapy. We recall that our 
model assumes the pressure to be homogeneous on the thorax. For each manipulation, we will focus on three quantity 
computed from the model outputs: 

• the quantity of mucus expelled out of the tree Vout - the net reduction of mucus volume in the lung after a 
manipulation. 

• the relative change in hydrodynamic resistance r during and after a manipulation - this quantity is correlated 
to the feeling for the patient to breath easily or not. 

• the mean mucus position in the tree, that reflects how deep mucus stands in the lung. It is equal to 


-”- + El'.o(^ v{si“-s+)+’) 


mmp = 
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All numerical experiments are performed with a background chest motion that corresponds to an idealized normal 
ventilation of the patient. We model normal ventilation as a sinusoidal negative pressure applied on the thorax with 
a frequency of 0.2 Hz and an amplitude = —5 cmH20, see figure [3]4: 

r,ventUr.: „ (1 - C0s(27rt/5)) 

^ext PI — „ 


With such data as an input, our model predicts a tidal volume of 0.5 L (6 L/min). This value for tidal volume is 
typical of rest regime ventilation [32] . 

The initial mucus distribution in the tree is plotted on figure [^. It is homogenous in the first six generations and 
fills about ten percents of the bronchi lumen area. Then mucus thickness increases in the three next generations 
to reach its maximal thickness in generation eight, filling about fifty percents of the bronchi lumen area. Finally, 
mucus thickness decreases and reaches zero percent in generation sixteen. The initial configuration for mucus has 
been chosen in such a way mucus is almost not affected by normal rest ventilation. 




generation index 


Figure 3. A: External oscillating pressure to model normal rest ventilation. B: Initial distribution of mucus in the tree: air 
is plotted in blue and mucus in green, bars height represents a proportion of airway surface at FRC. This distribution is not 
affected by normal rest ventilation. 


B. Manual physiotherapy 

Manual chest physiotherapy is performed by applying pressures on the chest with hands. The pressures are applied 
during the expiration phase only and practitioners are guided by chest motion and reaction to their manipulations. 
Although the hand locations play an important role in chest physiotherapy, our model accounts only for pressures 
amplitude and time dependence: pressure is applied homogeneously onto the chest. We mimicked chest physiotherapy 
in our model by altering the ventilation signal during expiration, as shown on figure |4]4: 

^ pventU^^^ ^ (sin(2^t/5), 0) 




Figure 4. A: Chest pressure t — > (t) used as an input for our model to mimic manual chest physiotherapy. Pressure 

is assumed homogeneous on the thorax, case with Pep = 20 cmH20. B: Relative hydrodynamic resistance r variations during 
a manipulation with Pep — 20 cmH20 (blue). The red line corresponds to FRC relative hydrodynamic resistance before 
the manipulation (r|t=os = 1.00), and the black line corresponds to FRC relative hydrodynamic resistance at the end of the 
manipulation (r|t =2303 = 0.79). 


Manipulation pressures are applied regularly during one session of 230 seconds, except for the 10 first and last 
seconds, where only normal ventilation occurs. The pressures applied are all the same during the whole manipulation. 
The twenty-five first seconds of pressure input are plotted on figure |4j4 with P^p = 20 cmH20. In manual chest 
physiotherapy, the pressure is not applied on the whole chest at once. The choice of the location for pressure 
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application is made by the physiotherapist in order to optimize the pressure effect. Our model is for now not able to 
account for pressure spatial inhomogeneities and as a first approximation, we assume the hand pressure spreads over 
the whole chest homogeneously. 

We investigated the role of pressure Pep amplitude on mucus motion and distribution and on hydrodynamic resis¬ 
tance changes. Smaller chest pressures induce smaller air velocities and thus smaller shear stresses in the bronchi. 
As a consequence, the stress threshold that has to be overcome for secretions to be mobilized is only reached in 
highly obstructed bronchi. But the motion is quickly stopped when a new equilibrium between air flow and secretions 
distribution is reached. In particular, for chest pressures lower than 16.5 cmH20, we do not observe in the model any 
secretions going out of the lung, see figure [5]4. This does not mean however that secretions distribution has not been 
affected, as shown on figure [^. An important point is that by performing the manipulation, secretions are moving 
in such a way that their distribution always decreases the total hydrodynamic resistance of the tree, down to a value 
corresponding to an equilibrium between secretions and air flows. The first steps of the manipulations are the more 
efficient, as shown on figure]^. The higher the pressure, the lowest is the tree hydrodynamic resistance at the end 
of the manipulation. The first conclusion is then that the manipulations decrease the hydrodynamic resistance of the 
patient, which may lead to an improvement of patient breathing quality, even if no secretions get out of the tree. 



Figure 5. A: Volume of secretions expelled at the end of the manipulation as a function of the pressure amplitude. B: Relative 
hydrodynamic resistance of the tree at the end of manipulation as a function of the pressure amplitude. C: Mean secretions 
position in the tree as a function of the pressure amplitude. Mean position is expressed in generation index. 


time: 0 s 


time: 25 s 






generation index generation index 


Figure 6. Evolution of mucus distribution in the bronchi during the manipulation, with Pep — 20 cmH20. The blue bars 
correspond to the surface of bronchi lumen area relatively to bronchi initial surface. The green bars correspond to the surface 
fraction of bronchi lumen area obstructed by secretions. At the end of manipulation, secretions have been globally motioned 
upward the tree, however some secretions have been pushed deeper in the lung, see generation 11 — 13. 


Secretions distribution in the tree are motioned from one generation to the next, and because manipulation pressures 
are applied during expiration, mucus tends to go up the tree, at least for moderate pressures amplitude, see figure 
[5p. The mechanism is as follow: during pressure application, secretions are moving upward in the tree from one 
generation to the next, and thus fill the upper generations. Lung’s recoil from the manipulation and lung opening 
due to inspiration are correlated to an inlet air flow in the lung, which potentially moves the secretions back into 
the lower generations. This backward secretions motion depends on how much secretions have accumulated during 
the pressure application. If the pressure is high, then secretions accumulate a lot in the upper bronchi, reducing 
drastically air lumen area, and by thus they are more affected by the inward air flow. This is why for high pressures, 
although there is less secretions in the tree at then end, they are globally deeper after the manipulation than before 
the manipulation. This phenomenon can be seen on figure where generations II to 13 have more secretions at the 
end of the manipulation. 
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C. Mechanical physiotherapy 

There exists a whole range of mechanical chest physiotherapy devices that are able to help secretions expectoration 
Elllllll]. We will focus on devices that work by applying oscillatory positive pressures on the chest, a method that 
is referred to as high frequency chest wall oscillation (HFCWO). HFCWO covers different ways to apply the pressure 
on the thorax. Chest compression (CC) devices [2] apply a global static pressure on the chest, and low amplitude 
oscillations are applied on the whole thorax, thus working at low lung volume. Focused pulses technique (FPT) devices 
[24] do not apply a global static pressure but instead apply high amplitude oscillations on specific locations on the 
chest using pistons. In this section, we propose to mimic with our model the effects of these two types of devices. We 
chose an input external pressure that adds to normal ventilation a static pressure Pg and an oscillatory pressure with 
amplitude Po and frequency /, see figure Q: 

= p:T“(t)+P- + Y “"( 2 ’'/*) 





Normal ventilation 


7X7^/V1 


B 


Figure 7. Chest pressure used as an input for our model to mimic mechanical chest physiotherapy based on high frequency 
chest wall oscillations (20 Hz in this example). Pressure is assumed homogeneous on the thorax. A: Pressure distribution on 
the thorax to mimic devices using chest compression (CC) using a static pressure. B: Pressure distribution on the thorax to 
mimic devices using focused pulses technique (FPT) without static pressure {Ps = 0). 


In the case of focused pulses technique, it is necessary to account for the pistons covering only a small fraction of the 
thorax. Thus, we use equivalent pressures Pq as inputs of the model. They represent the typical values applied by the 
pistons weighted by the surfaces ratio between pistons’ heads and thorax. As for manual chest physiotherapy modeling, 
manipulation pressures are applied for one session of 230 seconds. The 10 first and last seconds are normal ventilation. 
The time range 10s to 15s is used to start smoothly the manipulation by applying progressively the manipulation 
pressure. The time range 215s to 220s is used to stop smoothly the manipulation by decreasing progressively the 
manipulation pressure. 


1. An example of high frequency chest wall oscillation 

In this section, we study the predictions of the model with an input mimicking the effects of HFCWO with a static 
pressure Pg = 5.6 cmH20 and an oscillatory pressure Pq = 1.2 cmH20 at / = 20 Hz. Because of the static pressure, 
both ventilation and pressures oscillations are applied on a lung with low volume. Most particularly, this affects the 
net bronchi lumen area and the pleural pressure value and response mm- 

As for manual chest physiotherapy, our model predicts that the manipulation redistributes the secretions in the 
tree in such a way it reduces the hydrodynamic resistance of the tree by about twenty percents at the end of the 
manipulation, see figure [SjA. The manipulation efficiency is higher at the beginning and decreases with time, as shown 
on figure]^. Because the static pressure and the oscillations occurred during the whole ventilation cycle, the secretions 
spread a little more down the tree than for manual chest physiotherapy. Initially, mean mucus generation in the tree 
is 7.34, and it reaches 7.47 at the end of the manipulation, see figure]^ With the parameters used in this section, our 
model predicts that no secretions are going out of the tree. 


2. Role of static pressure Pg 

We investigated the role of the static pressure Pg that controls the lung volume at which the manipulation is 
performed. Oscillating pressure Pq was fixed to 1.2 cmH20 with a frequency / at 20 Hz. Surprisingly, we found 
that static pressure has only small effects on the final hydrodynamic resistance and on the mucus spreading in the 
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Figure 8. A: Evolution of mucus distribution in the bronchi during the manipulation with Ps = 5.6 cmH20 and Po = 
1.2 cmH20. B; Relative hydrodynamic resistance changes during a manipulation with Pep = 20 cmH20 (blue). The red line 
corresponds to FRC relative hydrodynamic resistance before the manipulation (rjt^os = 1.00) while the black line corresponds 
to FRC relative hydrodynamic resistance at the end of the manipulation (r|t= 230 s = 0.79). 


tree. Increasing the static pressure leads to a slight improvement of the hydrodynamic resistance at the end of the 
manipulation, see figure]^. Larger static pressure also sends secretions a bit deeper in the tree, and increases slightly 
the mean mucus generation at the end of the manipulation, see figure]^. Finally, increasing static pressure reduces 
to amount of secretions getting out of the tree, see figure [^. 



Figure 9. Role of static pressure Ps, with Po = 1.2 cmH20 and f — 20 Hz, data computed at the end of the manipulations for 
time t = 230 s. A: Volume of mucus getting out of the tree. B: Relative hydrodynamic resistance of the tree. C: Mean mucus 
position in the tree, in generation index. The red dashed line represents the value before the manipulation. 


These results indicate that, in the frame of our model, increasing the static pressure may have no real effect. Indeed, 
increasing the static pressure goes with a constriction of the bronchi, which increases the hydrodynamic resistance of 
the bronchial tree, as shown in the previous example on figure [^. The air flow created by the oscillating pressure 
is then smaller if the static pressure is larger. Shear forces determine the motion of the secretions, and in a bronchi, 
they are proportional to the ratio between the air flow rate in the bronchi and the diameter of the bronchi at the 
power three. Our results suggest that the increase on shear stress gained by smaller airways diameters is almost 
compensated by the loss in term of air flow rate. 


3. Influence of the oscillating pressure amplitude Po and frequency f 

Air flow in a bronchi determines the shear stress applied on the secretions and thus their motion. Air flow in 
a bronchi is partially related to velocity of the bronchi volume changes, see equation (|^. Thus secretions motion 
is correlated to bronchi volume change velocity. This quantity is actually the ratio between the amount of volume 
change over the time during which this volume change occurs. The amount of volume change is given by oscillating 
pressure amplitude, while the time for these changes to occur is given by the oscillating pressure frequency. Thus, 
these two quantities are ought to play an important role on manipulation efficiency. 

We investigated first the role of the oscillating pressure Po- Because static pressure has only a small effect in our 
model, we assume Pg = 0.6 cmH 20 which makes the total pressure applied on the thorax oscillate between 0 and 
1.2 cmH^O. The frequency / is kept at 20 Hz. As predicted, the larger the amplitude of the oscillations, the more 
efficient is the manipulation, see hgure |10[ As before, a more efficient mechanical manipulation tends to send more 
mucus down the tree, as shown on figure |10Tll, because oscillations are symmetric and applied during the whole cycle 
of ventilation. 

The frequency / also plays an important role on the manipulation efficiency. There is a minimal frequency under 
which the manipulations have low efficiency. Once this minimal frequency is reached, the hydrodynamic resistance 
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Figure 10. Role of oscillating pressure Po, with Pa = 0.6 cmH20 and / = 20 Hz, data computed at the end of the manipulations 
for time t = 230 s. A: Volume of mucus getting out of the tree. B: Relative hydrodynamic resistance of the tree. C: Mean 
mucus position in the tree, in generation index. The red dashed line represents the value before the manipulation. 


and the mean mucus position in the tree are only slightly affected by any supplementary increase of the frequency, as 
shown on figure [m Only the quantity of secretions that get out of the tree continue to increase regularly with the 
frequency, as shown on figure [Tl| 4, mainly because the secretions in the root of the tree (mimicking the trachea) are 
all the more mobilized than the frequency is high. 





Figure 11. Role of oscillating pressure frequency /, with Pa = 0.6 cmH20 and Po = 1.2 cmH20, data computed at the end 
of the manipulations for time t = 230 s. A: Volume of mucus getting out of the tree. B: Relative hydrodynamic resistance 
of the tree. C: Mean mucus position in the tree, in generation index. The red dashed line represents the value before the 
manipulation. 


D. Shrek number: a measure for chest physiotherapy efficiency 

The mucus behavior in a generation i can be predicted with the ratio between the air constraints on the mucus 
measured by the air shear stress Ti in an airway and the shear stress threshold tJo of the mucus. In an airway with 
radius and length li, the air shear stress is r* = ^/(7rrf). We define the instantaneous Shrek number Sh as 

the average of the ratios Ti/ao over the tree generations: 


Instantaneous Shrek ^5] number: iSh 


I 

N 


N-l 


E 




(16) 


The instantaneous Shrek number can be approximated with a formula easier to evaluate, using lung hydrodynamic 
resistance Raw and air flow at mouth level $ 0 , see details in appendix]^ This approximation is based on the hypothesis 
that the ratios of length over diameter of all the airways are equal to 3, which corresponds to the mean value measured 

m- 


iSh 


12Nao 


(17) 


The larger the instantaneous Shrek number, the more secretions are mobilized. An instantaneous Shrek number smaller 
than I does not mean that mucus is not mobilized, since Shrek number represents an average over the generations. 
It means however that in many generations, the mucus is probably not mobilized. However it is important to notice 
that airways impaired by mucus will have the highest Shrek numbers in the tree. Because chest physiotherapists 
are modulating the flow in their manipulation to keep mucus moving in the tree, they are indirectly and intuitively 
modulating the instantaneous Shrek number to keep it high enough to get a response from mucus. These modulations 
are paired with the real-time changes of lung hydrodynamic resistance. 

Finally, the instantaneous Shrek number can be reformulated as a function of the homogeneous external pressure 
used in our model. The new formulation involves quantities related to lung function and to manipulation: the 
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compliance of the lung Cl that relates the external pressure to the lung volume change, and the velocity with which 
changes in external pressure occurs APext/At: 


iSh. 


ChRaw ^Pex 
l2Nao At 


(18) 


This formula comes from the simple fact that air flow in the lung equals lung volume change AV^ divided by the 
duration for this change At. Also the compliance of the lung checks Cl = ^Vl/^P ext- Combining these equations 
leads to the alternative formulation for instantaneous Shrek number iSh in equation (18). 


To be able to compare different chest physiotherapy manipulations in the hypotheses of our model, we compute the 
Shrek number which is the time average of the instantaneous Shrek number. It defines a global mean efficiency for 
the manipulation, with the limit that any average may have. If the manipulation has a duration T, then: 


Sh = 7^ [ iSh{t)dt 

J Jo 


(19) 


The mean Shrek number associated to normal ventilation in a normal lung is about 0.2. 

It is important to highlight that Shrek number is only an indicator which aggregates many information, thus its use 
is correlated to an understanding of the physical phenomena involved in chest physiotherapy and to an understanding 
of its limitation due to the fact that it is an average. 
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Figure 12. Relative reduction of lung hydrodynamic resistance after 230 seconds manipulations (model predictions) as a function 
of the time averaged Shrek number. 

On figure |12[ we compare the efficiency in term of hydrodynamic resistance decrease for all the manipulations we 
studied in the previous sections. Most particularly, we can see that two manipulations with the same Shrek number 
have globally the same effect on hydrodynamic resistance. As expected, modulating the static pressure in chest 
compression techniques leads to very similar Shrek numbers, while manual and focused pulses techniques are able to 
cover a wide range of Shrek number by modulating their respective pressure amplitude. 


E. Comfort number: a measure for chest physiotherapy comfort 


The instantaneous comfort number of a manipulation is the lung tissue pressure shift relatively to normal ventilation. 
If at a given time of a manipulation, the volume of the lung is Vl while it would have been without the 

manipulation, then the instantaneous comfort number is: 


iCom = 


PussuejVL) - PussuejVr"'"') 

p {j/ventili 

-^tissuey^l^ ) 


( 20 ) 


The instantaneous comfort number is zero if only ventilation occurs. The higher the number, the larger are the 
mechanical constraints in the lung tissues. Notice that this number only measures the mechanical stress inside the 
lung and is not directly related to its hydrodynamic resistance. For example if lung tissues are suffering high negative 
pressure because of an high inflation, the associated lung hydrodynamic resistance is probably low but mechanical 
stress is high. 
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The comfort number of a manipulation with duration T is the mean value computed over the whole manipulation, 


1 F 

Com = — iCom(t)dt 
T Jo 


( 21 ) 


As a consequence, a strong pressure applied for a short time can lead to similar comfort number as a moderate or 
even small pressure applied for a long time. 


<> Chest physiotherapy (amplitude 0 to 3000 Pa) 

O FPT frequency (10 to 30 Hz) 

V FPT oscill. amplitude (10 to 100 Pa) 

x Chest compression static pressure (0 to 1000 Pa, 60 Pa oscill.) 



Figure 13. Comfort number computed from model predictions for the different chest physiotherapy techniques mimicked in this 
study. The higher the comfort number, the larger are the stress in the lung tissues. The comfort number for normal ventilation 
is zero. 

On figure |13[ the comfort number for different idealized chest physiotherapy techniques studied in this paper are 
plotted. Chest compression techniques with moderate static pressure are associated to high comfort numbers, mainly 
because the static pressure is applied constantly during the manipulation. Focused pulses technique are associated to 
low comfort numbers, because they are applying low pressures on chest. Notice that comfort number is only related 
to the net compressive forces induced into the lung by the manipulations, comfort number is not able to measure the 
role the frequency of pressure application may play on the patient comfort. 

We recall that in our model the pressure is applied onto the whole chest. Thus in the case of manual chest 
physiotherapy, comfort number has to be adapted to reflect the fact that the pressure is applied only on a subpart of 
the chest. In a first approximation, it can be done by weighting the pressure used to compute comfort number with 
the ratio between the surface on which the pressure is applied over the whole chest area. 


F. Model limitations 

Because the model described in this work contains the biophysical phenomena we identified as major in chest 
physiotherapy (’’minimal” model), we are confident that it is able to predict main behaviors of chest physiotherapy. 
However, our predictions can only be qualitative or comparative since the biophysical phenomena included are sim¬ 
plified versions of the ’’real” phenomena and are based on static mechanics. For example, secretions are modeled by 
a simple Bingham fluid, or secretions thickness is assumed homogeneous in a bronchus. Moreover, quantitative pre¬ 
dictions would need to include in the model some other phenomena which, although not strongly driving the system, 
might play a role on the quantitative data. Thus, any interpretation based on the results of our model has to be made 
in the frame of the model’s hypotheses. 

Our model is intrinsically non-linear. Non linearity in the model comes from two different aspects. First, in the way 
we deal with the mechanics of the lung and thorax. For this aspect, our model does not rely on a modeling process 
but instead on -non linear- data about lung mechanical inner pressure and about lung volumes (pressure curves from 
Agostoni et al. my Secondly, non linearity arises from the fluid mechanics of mucus which is based on Bingham fluid. 
Bingham fluid flows only if inner shear stress is above a threshold. Although non linear, secretion mechanics remains 
a basic simplification of real secretion mechanics, which behavior may be highly dependent on the local physical and 
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biological environment. Air fluid mechanics in our model is linear, and we ignored inertia and turbulence which may 
break the homogeneity of flows distribution in the tree, even if the tree has symmetric branching |21j . 

Asymmetric flow distribution can also be a consequence of an asymmetry in the geometry, either induced by 
asymmetric branching or by asymmetric secretions distribution in the tree. This last phenomenon may play a crucial 
role when a specific bronchus is more occluded by secretions than its neighbors: air flows might be redistributed towards 
the lowest resistive bronchial pathways and possibly miss the bronchi where secretions quantity is the highest. Chest 
physiotherapists are able to counteract this phenomena by applying pressures in a non homogeneous manner on the 
thorax. By assuming in our model that pressure is applied homogeneously on the thorax, none of these effects can be 
predicted with our current model. 


Chest physiotherapy is highly patient dependent: manipulation are based on real-time response to the praticien 
feeling (hand feeling, noises, mouth air flows, etc.). The data we used in this model for secretions properties [14] . 
pressure-volume relationships from Agostoni et al [T] and lung geometry and mechanics |16j reflect either mean or 
specific behaviors and cannot be applied as such to a specific patient. It has to be clear that the predictions of our 
model are general predictions that cannot be mapped as such to any patient. Again, careful interpretations in the 
frame of the model hypotheses have to be made before reaching any specific conclusion. 


We proposed two numbers, the Shrek number that measures an efficiency for the manipulation and the comfort 
number that measures a degree of compressive discomfort correlated to a manipulation. Both these numbers have 
been tested with our idealized model. They are however not validated in a medical way and there significance remains 
as of today an hypothesis. In this paper we propose these numbers as candidates for evaluating chest physiotherapy 
manipulations, but clinical measurements have to be made before any medical applications of these numbers are 
made. This leads to another question: how to evaluate quantitatively the values for these numbers in a clinical frame. 
Concerning Shrek number, we proposed a formula that allows to compute its value with typical medical data, see 
equation (18) and appendix]^ Concerning comfort number, it is necessary to get information about patient pressure- 
volume lung curves, which may be difficult in the frame of everyday consultations. This may be solved by assuming 
that patients pressure-volume lung curves have globally the same shapes and can be fitted through the measurement 
of a few simple physiological parameters at the beginning of a manipulation. This aspect is however out of the scope 
of this paper since it can be made only in a clinical frame. 


G. Conclusion 


We developed a model of the thorax and lung that is able to link the pressure on the thorax with the air and 
secretions motion in the airways. The goal of this model is to give indications on how chest physiotherapy may 
interact with the biomechanics of the thorax and of the lung. The model we built in this work is based on core 
phenomena that are involved in chest physiotherapy and can predict qualitative behaviors and compare different 
techniques, in the limitations of the model’s hypotheses. 

Our results indicate that manipulations need to overcome secretions threshold in order to be able to mobilize 
secretions. Our model shows that manipulations main effect is to reduce the hydrodynamic resistance of the airway 
tree by secretions redistribution in the tree, eventually reaching a distribution that is not anymore affected by the 
manipulation. This effect means that patient breathing might be instantly improved by the manipulation. The 
second main effect is mucus expectoration, which is also another mean to improve the status of the patient. Mucus 
expectoration predicted by our model is non negligible only for pressures high enough, and also, in the case of HFCWO, 
for frequencies high enough. 

All these phenomena are driven by the fact that the motor for secretions motion is the shear stress transmitted 
by air to the secretions. In a bronchus the shear stress is proportional to the air flow amplitude and to the inverse 
of the bronchus radius. Since air flow is created by applying pressure on the chest, the same method that is used 
to create the air flow will increase the hydrodynamic resistance of the airway tree, thus negatively affecting the air 
flow amplitude. Consequently, at some point, using a higher pressure will not really improve the efficiency of the 
manipulation, as shown on figure This phenomenon is well known by chest physiotherapists and is also at the 
origin of forced expiration curve shape |15| . 

Finally, we proposed two numbers to measure the efficiency of the manipulations: the Shrek number and the comfort 
number. They are both performing well in the case of our idealized manipulations. They need however to be validated 
with clinical studies. 
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Appendix A: Computation of pressure drop per unit length - function F 


In this appendix, we describe how we compute the pressure drop per unit length in a branch, C = dpjdz = 
F{^ai Sa, Sb), see equation (11). The way we compute the function F is identical to that used in [^. The expression 


of f is a consequence of air fluid dynamics and depends on the mucus states. 

All calculations take place in a branch of the tree whose radius is Sb/TT and air lumen area radius is 

TT. The yield stress of mucus is cto and its viscosity when it flows is pm- Mucus is solid on the range of 

radius [ra,ro] and liquid on the range [ro,rb]. If ro > rb then mucus is fully solid and if rp < Va then mucus is fully 
liquid. The yield radius ro is equal to |^|, see [5D]. Air viscosity is pa- 

We call C = dp/dz the pressure drop per unit length in the branch. The air flow in the branch is known and is 

equal to Finally, we call s the sign of C = dpjdz (s = 1 if C < 0 and s = —1 if C < 0). 

Because of developed flow and axi-symmetry hypotheses, fluid dynamics equations reduce to, see ^0] : 


1 d dp dv 

for r e [0,ra], (rS^r) + ^ = 0 with 

r dr dz dr 

By integration on r, the previous equation becomes: /iafy = 

The air flow in a branch is equal to: d>a = 27r v(r)rdr and, integrating by parts, it becomes: 


„2 4 


where v{ra) is the velocity of the air/mucus interface. This velocity depends on the state of mucus. 


First state: mucus is liquid between the branch wall located at Vb and the radius rg and solid elsewhere (case r^ < 
1^0 < fb)- In that case v{ra) = v(ro) since mucus is solid between rg and Tq. Then v(ra) = 'c(rg) = — 4 ^(ch—rg)^. 
Replacing v(ra) with the previous expression in the expression of and using the fact that rg = |2(Tg/C'|, then 
C is the solution of a second degree polynomial. The two solutions are: 


^ _ 1 STTrlpanao - 84>aPaPrn + -^TTrlplrbCJQC/apLm + mplpi/n " ^Tr'^r^PgalpLrn 

^ 2 + 2TTrlparl 

and 

^ ^ 1 STTrlpaTbao - ‘&4'aPaPm " 4x/-87rrg^^rbCro(^a/i^ + - 27r^rg^aggMm 

^ 2 + 2TTrlparl 

We have two possible pressure drops per unit length which gives two possible values for rg = |2CTg/C'|, however 
they are easily discriminated since one only at a time can check rj, < rg < r^. 

Second state: mucus is liquid everywhere, i.e. rg < r^. Then, v{ra) = — ^b)(?'a + ’"b ~ 2rg). Mixing with 

expression leads to 


rlprn + 2rlpa “ 2par'/ 

Third state: mucus is completely solid (rg > r^). In that case v{ra) = 0 and, 

27rr^ 


Th function F is built thanks to these formulas for C. To discriminate between the different possibilities and know 
the mucus state, we compute for each case the value(s) of C and compute its associated radius rg = |2(Tg/C'|. The 
case is correct only if the position of the radius rg relatively to r^ and rj, is compatible with the case hypotheses. 
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Appendix B: Computation of mucus flow - function G 


As for F, the computation of the mucus flow in a branch = G(C, Sa, Sb) depends on the state of mucus, see 
equation (14). The way we compute the function G is identical to that used in [20]. All calculations take place in a 


branch of the tree whose radius is Vb = and air lumen area radius is Va = \J S'o/tt. The yield stress of mucus is 

CTo and its viscosity when it flows is [im- Mucus is solid on the range of radius [va, ro] and liquid on the range [rp, rb\- 
If tq > Tb then mucus is fully solid and if Tq < then mucus is fully liquid. The yield radius vq is equal to I ^ I, see 


Fluid dynamics of mucus stands in the range [ra,rb] and reduces to 


^ = 0 where mucus is solid (typically in [ra,ro]) 

'Szr = ^ = o'o + where mucus is liquid (typically in [ro,r;i]) 


The flow of mucus is given by 


r” 

= 2 tt v{r)rdr 


Solving fluid dynamics equations for air and mucus then leads to analytical expressions for mucus flow in the branch: 

First case: mucus is liquid between the branch wall located at Vb and the radius tq and solid elsewhere (case < tq < 
rb). In that case v{ra) = v{ro) since mucus is solid between tq and r^. Then v{ra) = v{ro) = — 4 ^{ft — 
and 


= 27r^^^(r2 - r^) - 27r-^(rg - 

2 16 fi. 


G 


Second case: mucus is liquid everywhere, i.e. rg < Va- Then, v{ra) = 4 ^(ra — rb)(ra + rb — 2ro), and 

= -27r-^(rg -r2)2 


Third case: mucus is completely solid (rg > r^). In that case v(ra) = 0 and = 0. 

The value of rg relatively to Va and rb allows to find in which case we are and then to compute the function = 

G(G,Sa,Sb). 


Appendix C: Numerics 


The mathematical problem arising from the mechanics are solved using numerical computations. If n = (vi)i and 
w = (wi)i are two vectors then we define the product v * w as the vector (yiWi)i. If A is a matrix n x n and v a 
vector of length n then the matrix-vector product is written A.x. The equations are rewritten in a vectorial form 
using matrices-vector product . and vector-vector product 


Sb = H(Pext, Pair) 

^ Lb — '^alv ^alv 

C = F{^a.Sa.Sb) 


^ Pair — L.C 

^„,,o = G{G,Sa,Sb) 

^m.i = \ max (0, U.^rn,o) “ 2 min (0, H.d'm.o) 


Each vector X has its coordinates denoted XG) ^ for example Sb = (Sj^^ )z', the number of vector components is 
equal to the number of generations considered in the model (23). A, L, U and D are square matrices, see appendix 


|d for details on their expressions. G and F are applied to vectors in this way: G{G, Sa, Sb) = {G{GG 
and F{^a,Sa,Sb) = S^^\ si^'^))z. 


q { z ) q{z 

5 5 
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The equations lead to a non linear vectorial ordinary differential system. It is solved in a full implicit way. Given 
an initial condition at time 0 and the time-dependence of the tissue pressure t —>■ Pussueit), the system is discretized 
in this way: 

' Sbitn) = H{Pext{tn),Pair{tn)) 

$a(t„) = A-^ 

C{tn) = F{^a{tn),Sa{tn), Sb{tn)) 

^ Pairitn) = L.C{tn) 

Sb{tn) — Sa{tn) = Sb{tn-l) — S'a(tn-l) + (tn — tn-l)i\^m,l\ — |‘f’m,o|) * 

= G(C{tn), Sa{tn), Sb{tn)) 

. = 5 max( 0 , C/.$m,o(in))-2min(0,L».$m.o(i™)) 

Concretely, we are able to eliminate all -vectorial- variables except C{tn) and Sa(tn), thus this system reduces to 
two non linear vectorial equations: 

/ Sa{tn) = Pi {Sa(tn),C{tn), 5a (t„_ 1 ) , S';, (t„_ i ) ) 

\ C(t„) = P2 {Sa{tn),C{tn),Sa{tn-l)) 

The system is solved in this way at each time step n: 

• in a first step, we build a numerical function C„(5) such that C„(5) = P 2 {S,Cn{S), Saitn-i))- The function C„ 

is well defined since Cn'^ it is a strictly decreasing function of We use a Newton method to compute C„, 
initialized by Cinit = L.G(^,Sa{tn),S{Ptissue(tn) - Pair{tn-l))) with (j) = A-^ *Lby 

• in a second step, we solve the equation S = Pi {S,C{S), Sa{tn-i), Sb{tn-i)), again using a Newton method 
initialized with the air lumen area of the previous time step: Sinit = Sa{tn-i)- 

We use a parallel Newton method implemented in C-t--!- and OpenMP. All equations and variables are normalized. 


Appendix D: Computation of matrices for the problem vectorial formulation 


In this appendix, we give the expression of the matrices A, L, U and D used in the mathematical formulation of 
the model. In this appendix, the letter i and j refer to generation index and belong to the set {0,1, 2,..., N — 1}. 
The matrix L relates the pressure drop per unit length to the mean pressure inside the bronchi, see equation ( |12[ ): 

( 0 if j < f 
— \ if j — i 
y k A j > i 

The matrix U is used to compute the mucus flow when it goes up in the tree (toward the upper bronchi): 


= 


0 if j ^ i + 1 
1 if j = i -|- 1 


The matrix D is used to compute the mucus flow when it goes up in the tree (toward the upper bronchi): 


D 




0 if j ^ i — 1 
1 A j = i-l 


The matrix A relates the volume change of a bronchus with the air flows coming from the two daughter bronchi: 
A = Id- 2U. 


Appendix E: Approximation of the Shrek number 


In this appendix, we detail how the Shrek number can be approximated with equation ( |17[ ), starting from its 
definition in equation (16). First, we need to make the hydrodynamic resistance of an airway in generation i, 
Ri = 8pah/ appear in the sum: 


Sh = 


1 

N 


N-l 


E 


TrrfcTo 


4 

Nao 


N-l 


E 

2 = 0 


8k TTrf “ 
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Then, using the fact that the tree is symmetric, we have $ 


(i) 


Sh = 


1 

2Nao 


N-l 

E 


2 = 0 


= thus: 

k 2 * “ 


Because is the flow in the trachea, then it corresponds to the mouth flow 
g since it is the mean value measured in the lung [29) : 


Sh = 


12Nao 


N-l 


E 


Rj 


Finally, the sum ^ exactly the equivalent hydrodynamic resistance i?, 

symmetric branching [HI [22]. Finally: 


Sh = 


Raw^ 

12N<7 


a 

0 


The ratio ^ is approximated by 


of our model of the lung with 





